program eod_boundaries

    character :: clinoform_filename_P*1000,clinoform_filename_P1*1000, B*5, filepath*1000
    parameter (SPACE=2000)
    integer:: P, max,software_package, m, n, fb_number
    real :: xmin, xmax, ymin, ymax
    real :: EOD(2,SPACE,SPACE)
    real :: EOD_2(2,SPACE,SPACE)
    real :: EOD_3(2,SPACE,SPACE)
    real :: CLINOFORM_P(2,SPACE,SPACE)
    real :: CLINOFORM_P1(2,SPACE,SPACE)   
    
    write(*,*)"RMS press 1, Petrel press 2"
    read(*,*) software_package
    
    write(*,*)"How many facies boundaries? (max is 3)"
    read(*,*) fb_number
    
call facies_boundaries(EOD,EOD_2,EOD_3,software_package,fb_number, xmin, xmax, ymin, ymax)
call listing(SPACE,max)

    write(*,*)"Enter the file path for the clinoform surfaces (also the output destination)"
    read(*,*) filepath 

do P=1,(max-1)
    call clinoform_filename(clinoform_filename_P,clinoform_filename_P1,P, max, filepath)
    call clinoform_surfaces(CLINOFORM_P,CLINOFORM_P1,clinoform_filename_P,clinoform_filename_P1, n, m,P, xmin, xmax, ymin, ymax)
    call calculate_new_surfaces(CLINOFORM_P,CLINOFORM_P1,EOD,EOD_2,EOD_3,n,m,P,software_package,max, fb_number, xmin, xmax, ymin, ymax, filepath)
end do 

    open(unit=2, file="listing.txt")
    close(unit=2, STATUS='DELETE')
    open(unit=4, file="temp.txt")
    close(unit=4, STATUS='DELETE')
    open(unit=5, file="temp_2.txt")
    close(unit=5, STATUS='DELETE')     
    write(*,*) 'Press return to exit program'
    read(*,*)

stop
end program eod_boundaries
!----------------------------------------------------------------------------------------------
subroutine array_allocation (temp_name,TEMP,n,m, nonvalue, xmin, xmax, ymin, ymax)

    character :: temp_name*1000
    character :: dummy1*10, dummy2*10, dummy3*10, dummy4*10
    character :: dummy5*10, dummy6*10, dummy7*10, dummy8*10
    parameter (SPACE=2000)
    integer :: m, n, lines, rest, i, j, h, raw, mode, k
    real :: numb, nonvalue, xmin, xmax, ymin, ymax, depth_min, depth_max
    real :: passagefive(5),passagefour(4),passagethree(3),passagetwo(2)
    real :: passageone
    real :: passage(SPACE,5)
    real :: TEMP(2,SPACE,SPACE)
    
   open(unit=1, file=temp_name)
	 read(1,*) dummy1
	 read(1,*) dummy2
	 read(1,*) dummy3
	 read(1,*) dummy4
	 read(1,*) dummy5
	 read(1,*) dummy6
	 read(1,*) numb, nonvalue
	 read(1,*) m, n, xmin, xmax, ymin, ymax
	 read(1,*) dummy7
     read(1,*) dummy8
     
     do i=1, SPACE
       do j=1,SPACE
          TEMP(1,i,j)=nonvalue
          TEMP(2,i,j)=nonvalue
       end do
    end do
      
     lines=m/5+1
     rest=m-(lines-1)*5
	if (rest==0) then;
    	lines=lines-1
    end if
    
	 do i=1,n
    	h=1
        do j=1,(lines-1)
           read(1,*) passagefive
           passage(h,1:5)=passagefive
           h=h+1
    	end do

		if (rest==1) then;
          read (1,*) passageone
          passage(h,1)=passageone
        else if (rest==2) then;
          read (1,*) passagetwo
          passage(h,1:2)=passagetwo
        else if (rest==3) then;
          read (1,*) passagethree
          passage(h,1:3)=passagethree
        else if (rest==4) then;
          read (1,*) passagefour
          passage(h,1:4)=passagefour
        else;
          read (1,*) passagefive
          passage(h,1:5)=passagefive
		end if
        
		do k=1,m
          raw=k/5
          mode=k-raw*5
          if (mode==0) then;
            mode=5
            raw=raw-1
          end if
          TEMP(1,i,k)=passage(raw+1,mode)          
         end do
     end do
    close(unit=1)   
              
return
end subroutine array_allocation
!---------------------------------------------------------------------------------------------------
subroutine facies_boundaries (EOD,EOD_2,EOD_3,software_package,fb_number, xmin, xmax, ymin, ymax)

    character :: temp_name*1000
    parameter (SPACE=2000)
    integer :: n, m, i, j, k,software_package, fb_number
    real :: nonvalue, xmin, xmax, ymin, ymax
    real :: TEMP(2,SPACE,SPACE)
    real :: EOD(2,SPACE,SPACE)
    real :: EOD_2(2,SPACE,SPACE)
    real :: EOD_3(2,SPACE,SPACE)
  
             
    write(*,*)"Enter the file name of the base of EOD 1:"
    read(*,*) temp_name
    call array_allocation (temp_name,TEMP,n,m, nonvalue, xmin, xmax, ymin, ymax)
     do i=1,n
        do j=1,m
           EOD(1,i,j)=TEMP(1,i,j) 
        end do
     end do
     !k=1
     !call write_zmap_file(EOD,temp_name,n,m,k, software_package, xmin, xmax, ymin, ymax)  
     
     
  if (fb_number >1) then;
     
    write(*,*)"Enter the file name of the base of EOD 2:"
    read(*,*) temp_name
    call array_allocation (temp_name,TEMP,n,m, nonvalue, xmin, xmax, ymin, ymax)
     do i=1,n
        do j=1,m
           EOD_2(1,i,j)=TEMP(1,i,j) 
        end do
     end do
     !k=1
     !call write_zmap_file(EOD_2,temp_name,n,m,k, software_package, xmin, xmax, ymin, ymax)
     
   end if
 
  if (fb_number >2) then;
    
    write(*,*)"Enter the file name of the base of EOD 3:"
    read(*,*) temp_name
    call array_allocation (temp_name,TEMP,n,m, nonvalue, xmin, xmax, ymin, ymax)
     do i=1,n
        do j=1,m
           EOD_3(1,i,j)=TEMP(1,i,j) 
        end do
     end do
     !k=1
     !call write_zmap_file(EOD_3,temp_name,n,m,k, software_package, xmin, xmax, ymin, ymax)
      
   end if

return
end subroutine facies_boundaries
!----------------------------------------------------------------------------------------------------
subroutine listing(SPACE,max)

 integer :: i, max

 open(unit=2, file="listing.txt")
  do i=1,SPACE
   write(2,*) i
  end do
 close(unit=2)
 
 write(*,*)"How many clinoform surfaces are there?"
 read(*,*) max

 return
end subroutine listing
!----------------------------------------------------------------------------------------------------
subroutine clinoform_filename(clinoform_filename_P,clinoform_filename_P1, P, max, filepath)

    character :: clinoform_filename_P*1000,clinoform_filename_P1*1000,B*5, temp_name*1000, filepath*1000 
    integer :: i, P, max
          
    open (unit=2,file="listing.txt")
        do i=1,P
	       read(2,*) B
        end do
    close(unit=2) 
    clinoform_filename_P=trim(filepath)//"Clinoform_Surface"//trim(B)//".zmap+"
  

    open (unit=2,file="listing.txt")
        do i=1,P+1
	       read(2,*) B
        end do
    close(unit=2)  
    clinoform_filename_P1=trim(filepath)//"Clinoform_Surface"//trim(B)//".zmap+"
      
return
end subroutine clinoform_filename
!--------------------------------------------------------------------------------------------------------
subroutine clinoform_surfaces (CLINOFORM_P, CLINOFORM_P1,clinoform_filename_P,clinoform_filename_P1, n, m, P, xmin, xmax, ymin, ymax)

    character :: clinoform_filename_P*1000,clinoform_filename_P1*1000,temp_name*1000, B*5 
    parameter (SPACE=2000)
    integer :: n, m, i, j, k, P, max, software_package 
    real :: nonvalue, xmin, xmax, ymin, ymax
    real :: TEMP(2,SPACE,SPACE)
    real :: CLINOFORM_P(2,SPACE,SPACE)
    real :: CLINOFORM_P1(2,SPACE,SPACE)
    
    temp_name=clinoform_filename_P
    call array_allocation (temp_name,TEMP,n,m, nonvalue, xmin, xmax, ymin, ymax)
     do i=1,n
        do j=1,m
           CLINOFORM_P(1,i,j)=TEMP(1,i,j)
        end do
     end do
     k=1
         
    temp_name=clinoform_filename_P1
    call array_allocation (temp_name,TEMP,n,m, nonvalue, xmin, xmax, ymin, ymax)
     do i=1,n
        do j=1,m
           CLINOFORM_P1(1,i,j)=TEMP(1,i,j)
        end do
     end do
     k=1
        
return
end subroutine clinoform_surfaces
!----------------------------------------------------------------------------------------------------
subroutine write_zmap_file(TEMP,temp_name,n,m,k, software_package, xmin, xmax, ymin, ymax)

    character :: temp_name*1000
    character :: path*97
    parameter (SPACE=2000)
    integer :: n, m, i, k, software_package
    real :: TEMP(2,SPACE,SPACE)
       
    !!path="C:\PhD\Algorithms\Programs\Clinoform Algorithm_V6\clinoform algorithm v6\clinoform algorithm v6\PSS1 Surfaces\PS1.5\Facies Boundaries"

     open(unit=3,file=temp_name)
        write(3,*)'!'
        write(3,*)'!     ZIMS FILE NAME :  Envtopsphere.txt'
        write(3,*)'!     FORMATTED FILE CREATION DATE: ??/??/2010'
        write(3,*)'!     FORMATTED FILE CREATION TIME: ??:??'
        write(3,*)'!'
        write(3,1)'@Grid HEADER, GRID, 5'
1       format(A21)
        write(3,*)'	15,-9999.000,  ,	4 ,1.0000'
        write(3,"(I6,A,I5,A,1F12.4,A,1X,1F12.4,A,1F12.4,A,1F12.4)") m,',',n,',',xmin,',',xmax,',',ymin,',',ymax
        write(3,*)'0.0000     0.0000     0.0000'
	    write(3,2)'@'
2       format(A1)
        do i = 1, n 
            write(3,"(1F10.4,1X,1F10.4,1X,1F10.4,1X,1F10.4,1X,1F10.4)") TEMP(k,i,1:m)
         end do
    close(unit=3)

return
end subroutine write_zmap_file
!--------------------------------------------------------------------------------------------------------------------------------
subroutine calculate_new_surfaces(CLINOFORM_P,CLINOFORM_P1,EOD, EOD_2,EOD_3,n, m, P, software_package,max, fb_number, xmin, xmax, ymin, ymax, filepath)

    character :: temp_name*1000, G*5, W*5, filepath*1000
    parameter (SPACE=2000)
    integer :: P, i, j, n, m, k, software_package, u, max, c, fb_number, q
    real :: xmin, xmax, ymin, ymax
    real :: EOD(2,SPACE,SPACE)
    real :: EOD_2(2,SPACE,SPACE)
    real :: EOD_3(2,SPACE,SPACE)
    real :: CLINOFORM_P(2,SPACE,SPACE)
    real :: CLINOFORM_P1(2,SPACE,SPACE)
    real :: INT_GRIDDING_SURFACE(2,SPACE,SPACE)
    real :: INT_GRIDDING_SURFACE_2(2,SPACE,SPACE)
    real :: INT_GRIDDING_SURFACE_3(2,SPACE,SPACE)
    
    u=(max-P)
    open(unit=4,file="temp.txt", STATUS='REPLACE')        
	      write(4,*) u
    close(unit=4)
    
    open(unit=4,file="temp.txt") 
    	 do i=1,1
	      read(4,*) G
	    end do
    close(unit=4)
    
    do i=1,n
      do j=1,m
         INT_GRIDDING_SURFACE(1,i,j)=-9999.0000
         INT_GRIDDING_SURFACE_2(1,i,j)=-9999.0000
         INT_GRIDDING_SURFACE_3(1,i,j)=-9999.0000
      end do
    end do
    
    do i=1,n
       do j=1,m
          
          if(software_package==1.and.CLINOFORM_P1(1,i,j)>=EOD(1,i,j))then;
             INT_GRIDDING_SURFACE(1,i,j)=CLINOFORM_P1(1,i,j)
          else if (software_package==2.and.CLINOFORM_P1(1,i,j)<=EOD(1,i,j))then;
            INT_GRIDDING_SURFACE(1,i,j)=CLINOFORM_P1(1,i,j)
          end if
          if(software_package==1.and.CLINOFORM_P1(1,i,j)<=EOD(1,i,j))then;
              INT_GRIDDING_SURFACE(1,i,j)=EOD(1,i,j)
          else if(software_package==2.and.CLINOFORM_P1(1,i,j)>=EOD(1,i,j))then;
             INT_GRIDDING_SURFACE(1,i,j)=EOD(1,i,j)
          end if
          if(software_package==1.and.CLINOFORM_P(1,i,j)<=EOD(1,i,j))then;
             INT_GRIDDING_SURFACE(1,i,j)=CLINOFORM_P(1,i,j)
          else if(software_package==2.and.CLINOFORM_P(1,i,j)>=EOD(1,i,j))then;
             INT_GRIDDING_SURFACE(1,i,j)=CLINOFORM_P(1,i,j)
          end if
   
       end do
    end do
    
   
      k=1
      temp_name=trim(filepath)//"Clinoform_Surface"//trim(G)//".FB1.zmap+"
   call write_zmap_file(INT_GRIDDING_SURFACE,temp_name,n,m,k, software_package, xmin, xmax, ymin, ymax)     

   
if (fb_number>1) then;

   do i=1,n
       do j=1,m
          
          if(software_package==1.and.CLINOFORM_P1(1,i,j)>=EOD_2(1,i,j))then;
             INT_GRIDDING_SURFACE_2(1,i,j)=CLINOFORM_P1(1,i,j)
          else if(software_package==2.and.CLINOFORM_P1(1,i,j)<=EOD_2(1,i,j))then;
             INT_GRIDDING_SURFACE_2(1,i,j)=CLINOFORM_P1(1,i,j)
          end if
          if(software_package==1.and.CLINOFORM_P1(1,i,j)<=EOD_2(1,i,j))then;
              INT_GRIDDING_SURFACE_2(1,i,j)=EOD_2(1,i,j)
          else if(software_package==2.and.CLINOFORM_P1(1,i,j)>=EOD_2(1,i,j))then;
                INT_GRIDDING_SURFACE_2(1,i,j)=EOD_2(1,i,j)
          end if
          if(software_package==1.and.CLINOFORM_P(1,i,j)<=EOD_2(1,i,j))then;
             INT_GRIDDING_SURFACE_2(1,i,j)=CLINOFORM_P(1,i,j)
          else if(software_package==2.and.CLINOFORM_P(1,i,j)>=EOD_2(1,i,j))then;
              INT_GRIDDING_SURFACE_2(1,i,j)=CLINOFORM_P(1,i,j)
          end if
   
       end do
    end do
    

      k=1
      temp_name=trim(filepath)//"Clinoform_Surface"//trim(G)//".FB2.zmap+"
   call write_zmap_file(INT_GRIDDING_SURFACE_2,temp_name,n,m,k, software_package, xmin, xmax, ymin, ymax)     

   
   do i=1,n
      do j=1,m
          
          if(software_package==1.and.CLINOFORM_P1(1,i,j)>=EOD_3(1,i,j))then;
             INT_GRIDDING_SURFACE_3(1,i,j)=CLINOFORM_P1(1,i,j)
          else if(software_package==2.and.CLINOFORM_P1(1,i,j)<=EOD_3(1,i,j))then;
             INT_GRIDDING_SURFACE_3(1,i,j)=CLINOFORM_P1(1,i,j)
          end if
          if(software_package==1.and.CLINOFORM_P1(1,i,j)<=EOD_3(1,i,j))then;
              INT_GRIDDING_SURFACE_3(1,i,j)=EOD_3(1,i,j)
          else if(software_package==2.and.CLINOFORM_P1(1,i,j)>=EOD_3(1,i,j))then;
              INT_GRIDDING_SURFACE_3(1,i,j)=EOD_3(1,i,j)
          end if
          if(software_package==1.and.CLINOFORM_P(1,i,j)<=EOD_3(1,i,j))then;
             INT_GRIDDING_SURFACE_3(1,i,j)=CLINOFORM_P(1,i,j)
          else if(software_package==2.and.CLINOFORM_P(1,i,j)>=EOD_3(1,i,j))then;
             INT_GRIDDING_SURFACE_3(1,i,j)=CLINOFORM_P(1,i,j)
          end if
   
       end do
    end do
    

      k=1
      temp_name=trim(filepath)//"Clinoform_Surface"//trim(G)//".FB3.zmap+"
   call write_zmap_file(INT_GRIDDING_SURFACE_3,temp_name,n,m,k, software_package, xmin, xmax, ymin, ymax)     

 !rename clinoform surfaces in correct startigraphic order for RMS gridding  
        !q=(max-P)+1
            !open(unit=5,file="temp_2.txt", STATUS='REPLACE')        
	            !write(5,*) q
            !close(unit=5)
    
            !open(unit=5,file="temp_2.txt") 
    	        !do i=1,1
	                !read(5,*) W
	            !end do
            !close(unit=5)
   

            k=1
            temp_name=trim(filepath)//"Clinoform_Surface"//trim(G)//".Clinoform.zmap+"
         call write_zmap_file(CLINOFORM_P1,temp_name,n,m,k, software_package, xmin, xmax, ymin, ymax)     

   
end if
     
   print*, "Completed clinothem",P
            
return
end subroutine calculate_new_surfaces